\begin{verbatim}
function RelativeDifferentialProbabilities(D1,D2,PT)                          
    // Given 2-character strings D1 and D2 representing
    // XY-ZW and XY-QR, return the 26x26 matrix of 
    // probabilities, (P(XY)P(ZW)P(QR)), scaled to a 
    // probability function on the set {A,..,Z}^2.
    // PT is a sample plaintext for use in determining
    // the 2-character frequency distribution of English.
    assert #D1 eq 2 and #D2 eq 2;                                             
    AZ := {@ CodeToString(64+i) : i in [1..26] @};                            
    r1 := Index(AZ,D1[1]); s1 := Index(AZ,D1[2]);                             
    r2 := Index(AZ,D2[1]); s2 := Index(AZ,D2[2]);                             
    FDD := RealField()!0;                                                     
    DD2 := MatrixAlgebra(RealField(),26)!0;                                   
    F2D := DigraphFrequencyDistribution(PT);                                  
    for i1, j1 in [1..26] do                                                  
        i2 := ((i1-r1) mod 26) + 1;                                           
        j2 := ((j1-s1) mod 26) + 1;                                           
        i3 := ((i1-r2) mod 26) + 1;                                           
        j3 := ((j1-s2) mod 26) + 1;                                           
        F3 := F2D[i1,j1] * F2D[i2,j2] * F2D[i3,j3];                           
        DD2[i1,j1] +:= F3;                                                    
        FDD +:= F3;                                                           
    end for;                                                                  
    return (1/FDD)*DD2;                                                       
end function;                                                                 
\end{verbatim}
